#include "ClassMC.hpp"

//  This funtion will be called inside likelihood_ClassMC(**)

double Likelihood_LSST_WB(  imcmc_double&   full_params,
                            double&         lndet,
                            double&         chisq,
                            void*           model,
                            void*           data,
                            istate&         state ){
    lndet = chisq = 0;

    LSST_Forecast_WB *lsst = static_cast<LSST_Forecast_WB*>(data);

    arma::vec Params = arma::zeros(11);
    // arma::vec Params = arma::zeros(10);
    Params(0) = full_params["w0_fld"];
    Params(1) = full_params["wa_fld"];

    // Omega_m   = full_params["Omega_m"]; // Omega_m is derived parameter.
    // H0        = full_params["H0"];
    // Params(2) = Omega_m*(H0/100.)*(H0/100.);
    Params(2) = full_params["omegamh2"];
    Params(3) = full_params["omegabh2"];
    Params(4) = full_params["theta_s"];
    Params(5) = full_params["omegak"];  // be careful!!
    Params(6) = full_params["tau_reio"];
    Params(7) = full_params["YHe"];
    Params(8) = full_params["n_s"];
    Params(9) = full_params["alpha_s"];
    // Params(10)= full_params["ln(10^10A_s)"]-_10_log10_;
    Params(10) = full_params["lnAs"];

    arma::vec dP = Params - lsst->Params_fid;
    chisq = arma::as_scalar( dP.t() * lsst->FisherMatrix * dP );

// Check parameter updating status. This part should be kept.
// cout << "debug:\n";
// cout << "Params_fid =\n" << lsst->Params_fid.t() << "\n";
// cout << "Params =\n" << Params.t() << "\n";
// cout << "dP =\n" << dP.t() << "\n";
// cout << " -- CHISQ = " << chisq << endl;

//  TEST the chisq calculation:
    // double chisq_temp = 0.0;
    // double *x = new double[11];
    // double *y = new double[11];
    //
    // for( int i=0; i<11; ++i ){
    //     x[i] = dP[i];
    // }
    //
    // for( int i=0; i<11; ++i ){
    //     y[i] = 0.0;
    //     for( int j=0; j<11; ++j ){
    //         y[i] += x[j]*lsst->FisherMatrix(j,i);
    //     }
    // }
    //
    // for( int i=0; i<11; ++i ){
    //     chisq_temp += x[i] * y[i];
    // }
    //
    // delete[] x;
    // delete[] y;

    // cout << "@debug chisq_arma = " << chisq << endl;
    // cout << "@debug chisq_test = " << chisq_temp << endl;
    // cout << "@debug chisq_diff = " << chisq_temp - chisq << endl;

    return -0.5*chisq;
}
